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We study the limiting behavior of large-amplitude standing waves on deep water using high- 
resolution numerical simulations in double and quadruple precision. While periodic traveling waves 
approach Stokes's sharply crested extreme wave in an asymptotically self-similar manner, we find 
that standing waves behave differently. Instead of sharpening to a corner or cusp as previously 
conjectured, the crest tip develops a variety of oscillatory structures. This causes the bifurcation 
curve that parametrizes these waves to fragment into disjoint branches corresponding to the different 
oscillation patterns that occur. In many cases, a vertical jet of fluid pushes these structures upward, 
leading to wave profiles commonly seen in wave tank experiments. Thus, we observe a rich array of 
dynamic behavior at small length scales in a regime previously thought to be self-similar. 



Singularities in fluid mechanics are generally expected 
to be asymptotically self-similar lj. These can be dy- 
namic singularities, such as bubble pinch-off [2] or wave 
breaking [3], or parametric singularities, where a family 
of smooth solutions terminates at a singular solution. A 
famous example of the latter type was posed by Stokes 
in 1880, who used an asymptotic expansion of the stream 
function to argue that the periodic traveling water wave 
of greatest height should have an interior crest angle of 
120°. This crest angle has been confirmed in numerous 
computational studies [I] as well as theoretically [5] . The 
asymptotic behavior of the almost highest traveling wave 
was analyzed by Longuet-Higgins and Fox [6] [7] . 

Because genuine dynamics are involved, existing nu- 
merical methods have been unable to maintain the ac- 
curacy needed to fully explore the limiting behavior of 
large-amplitude standing waves. As a result, Penney and 
Price's conjecture [8] that a limiting standing wave ex- 
ists and develops 90° interior crest angles each time the 
fluid comes to rest has remained open since 1952. Such 
a singularity would be both dynamic and parametric. 
The standing waves in question are spatially periodic and 
have zero impulse (horizontal momentum), maintaining 
even symmetry for all time. They are also temporally pe- 
riodic, alternately passing through two zero- velocity rest 
states of maximal potential energy. 

Small-amplitude standing waves of this type were 
proved to exist by Iooss, Plotnikov, and Toland [S]- 
Larger-amplitude waves were computed by Mercer and 
Roberts [10], who discovered that the wave steepness 
(half the crest-to-trough height) does not increase mono- 
tonically over the entire one-parameter family of standing 
waves. They proposed using (downward) crest accelera- 
tion, A c , as a continuation parameter instead. We re- 
produce (and extend) their plot of wave steepness versus 
crest acceleration in Fig. [T] Since pressure increases with 
depth near the free surface |llj . Euler's equations imply 
that A c cannot exceed g, the acceleration of gravity. 

Taylor |12j performed wave tank experiments and con- 
firmed that large-amplitude standing waves do form rea- 
sonably sharp crests close to 90 degrees. A further in- 
crease in amplitude caused the waves to splash and be- 



come unstable in the transverse direction. Grant [13] and 
Okamura [14] have written theoretical papers to support 
the 90° conjecture. Okamura also performed numerical 
experiments [151 116) to back this claim. Extrapolating 
from numerical solutions, Mercer and Roberts |10j spec- 
ulated that the limiting crest angle might be as sharp as 
60°. Schultz et. al. [17] also predicted a limiting wave 
profile with a crest angle smaller than 90° and offered 
the possibility that a cusp may form instead of a corner. 

Our objective is to challenge the assumption that 
standing waves behave as traveling waves in their ap- 
proach of an "extreme" limiting wave. If there is no 
limiting wave profile, then a local analysis suggesting a 
geometric singularity (corner or cusp) is inapplicable. 

The equations of motion for a two-dimensional irrota- 
tional ideal fluid of infinite depth are 

r)t = <Py - Vx<Px, (la) 
$ t =P[<f> y r, t -l\V<f>\ 2 -gr,], (lb) 

where rj{x, t) is the upper boundary of the evolving fluid 
and &(x, t) = (f>(x, rj(x, t), t) is the restriction of the veloc- 
ity potential to the free surface. Both rj(x,t) and &(x,t) 
are assumed to be 27r periodic in a;. In (flb|, P is the or- 




FIG. 1. Bifurcation diagram and selected standing waves, 
plotted at equal time slices over a quarter period. The wave- 
length is taken to be 2tt, and g = 1. The crest tip sharpens as 
Ac increases over the range < A c < 0.985, where previous 
numerical studies are reliable. In particular, the curvature at 
the crest is visibly higher for solution B than for A. 
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thogonal projection to zero mean. This equation comes 
from $ t = (j)t + 4>yT)t and the unsteady Bernoulli equation 
4>t + \ |"V0| 2 + - p +gy = c(t), where the arbitrary constant 
c(t) is chosen to preserve the mean of $(x, t). 

To evaluate the right-hand side of for the purpose 
of time stepping, we use a boundary integral collocation 
method. Details will be given elsewhere [18] . Briefly, we 
represent ^ at a point z — x + iy in the fluid using a 
double layer potential. Suppressing t in the notation and 
summing over periodic images |19j . the result is 
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K(z, a)u(a) da, 
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where K(z, a) — Im { ^ cot ( z 2 ) }• prime repre- 
sents a derivative with respect to a, and 



(3) 



is a parametrization of the curve. The change of variables 
x = £,(a) allows for smooth mesh refinement near the 
crest tip. Letting z approach the boundary, we obtain a 
second-kind Fredholm integral equation for /1: 
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K{a,P)n{P)dp, (4) 
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Once n{a) is known, we compute 4> x and 0^ on the 
boundary from closing the system ([I]); see [T8] . 

We discretize space and time adaptively to resolve the 
solution as it becomes increasingly singular. Time is di- 
vided into v segments 6{T, where 9\ + • • • + 6 V = 1/4 and 
T is the current guess for the period. On segment I, we 
fix the number of (uniform) time steps, Ni, the number 
of spatial grid points, Mi, and the function 



6(a) = 



/ Et{p)dp, E l {a) = l-P[A l sin 4 (a/2)] 
Jo 



which controls the grid spacing in the change of variables 
x = £i(a). Ai is a parameter chosen between (uniform 
spacing) and 8/5, the value where &(a) ceases to be a 
diffcomorphism. As before, P projects out the mean. 

To compute standing waves, we use the Levenberg- 
Marquardt method 2Q , a trust-region algorithm for non- 
linear least squares problems, to minimize 



f( c ) = J 7r f o Hx,T/4) 2 dx, cG 
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where c contains the period as well as the nonzero Fourier 
modes of the initial conditions; i.e., T = Cq and 

Vk(0) = c\ k \ (k odd), $ fc (0) = c\ k \ (k even). (6) 

Here k ranges from —n to n, excluding 0, and n is chosen 
to be close to \M\, leaving the upper half of the spectrum 



of 77 and $ to be zero initially. A symmetry argument 
[TO] shows that driving the velocity potential to zero at 
time T/4 with initial conditions of the form ^ leads to 
a standing wave with period T and zero impulse. The 
method fails if / reaches a nonzero local minimum. 
We discretize ^ with spectral accuracy by redefining 



where r £ E m , m = M„, and 



n = $(£ v (ai),T/4)y/E v (ai)/m, a t = 2ni/M v . 

The square root comes from dx — E v (a) da. Typically, 
An < m < lOn. To track families of solutions, one of 
the Cfe is chosen as a continuation parameter |21) and 
eliminated from the search space when minimizing /. 
When a turning point is detected in this c&, we switch 
to a different one; see [TBI US] f° r details. The Jacobian 
Jik = dri/dck is computed by solving the linearization 
of (JIJ) about the current solution to obtain g^-$(x,T/4). 
This can be parallelized very efficiently |18| , dramatically 
increasing the resolution we are able to achieve. 

Our results are summarized in Figs. [2] and [3] First, 
we corroborate the result of Mercer and Roberts [10] 
that wave steepness, h, reaches a local maximum of 
hmax = 0.62017 at A c = 0.92631. (The values reported in 
[TO] were 0.6202 and 0.9264.) Using quadruple precision, 
we are able to compute h max to 26 digits of accuracy 
and the corresponding A c to 13 digits. Okamura [15j . 
who found that h increases monotonically all the way to 
A c = 1, was incorrect. Second, we find that crest accel- 
eration has turning points at A c = 0.99135 and 0.99040. 
This is a surprise, as A c was chosen as a continuation 
parameter in [TO] to avoid the lack of monotonicity in 
h. In our work, h and A c are plotted parametrically as 
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FIG. 2. The bifurcation curve in Fig. [TJ becomes fragmented 
in the range 0.985 < A c < 1, where previous numerical stud- 
ies break down. The labels A-0 correspond to wave profiles 
shown in Fig. [3] The turning point in wave steepness at C, the 
lack of monotonicity in A c , the complicated branching struc- 
ture, and the existence of standing waves with h > 0.62017 
were not previously known. 
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FIG. 3. Evolution of standing waves and velocity potential over a quarter period. (Top left) Solutions A-O in Fig. [2] are 
plotted on top of each other at the indicated times. (Right) These solutions develop oscillatory structures near the crest that 
change phase across disconnections in the bifurcation diagram. Solutions D-O have 350-600 grid points between 0.99-7T and 
1.0l7r. With at most 3 grid points in this interval, previous numerical studies could not resolve these structures. (Bottom left) 
The velocity potential of solution O has been driven almost to zero at t = T/4, yielding / = 1.3 x 10 -26 . For this solution, we 
used v = 4, 6i = {0.2, 0.2, 0.4, 0.2}, M t = {4608, 6144, 6912, 8192}, iV ; = {192, 288, 768, 480}, and A\ = {0, 0.774, 1.358, 1.381}. 



functions of whichever Ck is currently used as a continu- 
ation parameter. Finally, in the process of tracking this 
primary branch of solutions, we discovered several other 
families of standing waves. Each of these branches was 
tracked in both directions until the computations became 
too expensive to continue further with the desired accu- 
racy, / ~ 10 -26 in double precision. 

The standing waves that constitute these branches look 
qualitatively similar to each other in the large, where they 
closely resemble the photographs from Taylor's wave tank 
experiments [T^]. However, as illustrated in Fig. [3j so- 
lutions on different branches feature different oscillation 
patterns in the vicinity of the crest tip. The rapid in- 
crease in wave steepness from solution E to solution O in 
Fig. [^corresponds to a vertical jet of fluid that forms near 
the crest before the standing wave reaches its rest state. 
The resulting protrusion causes the maximum slope to 
be much larger than 1 for most of these solutions. Taylor 
photographed similar structures at the crest in his wave 
tank experiments. Schultz et. al. [17] argued that surface 
tension was responsible for these protrusions, but we find 
that they occur even without surface tension. Compar- 
ing solutions A-E on the primary branch, we see that 
solutions eventually flatten out at the crest and become 
oscillatory rather than sharp. Figure [4] provides further 
evidence that these oscillations grow large enough to pre- 
vent this family of solutions from approaching a limiting 
wave profile in an asymptotically self-similar fashion. 

Regarding accuracy, our method is spectrally accurate 
in space, 8th or 15th order in time [TH], and quadratically 
convergent in the search for a minimizer of / in ([5|. We 
achieve robustness by formulating the shooting method 
as an overdetermined nonlinear least squares problem. 
If the numerical solution loses resolution, the equations 
7"i(c) = become incompatible with each other and the 



objective function / = \r T r grows accordingly. This pre- 
vents the method from giving misleading overestimates of 
the accuracy of the standing waves it finds. For example, 
we recomputed solution O of Fig. [3] in quadruple preci- 
sion on a finer mesh (M x = {6144, 7500, 8192, 9216}, 9 t = 
{0.1,0.3,0.4,0.2}, and A x = {0,1.043,1.405,1.476}), us- 
ing the initial conditions obtained by minimizing / in 
double precision. The more accurately computed value 
of / is 8.6 x 10~ 27 , which is 34% smaller than predicted 
in double precision. This level of inaccuracy in the pre- 
dicted error is acceptable, as driving T/4) to zero 
entails eliminating as many significant digits as possible. 
For solution A, we repeated the minimization in quadru- 
ple precision, causing / to decrease from 2.2 x 10 -28 to 
2.1 x 10~ 60 . In addition to /, we monitor energy conser- 
vation and the decay of Fourier modes at various times 
to ensure that r] and <£> remain resolved to machine pre- 
cision; see [H] for more details. 
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FIG. 4. Breakdown of self-similarity on the primary branch. 
When lengths are rescaled so the radius of curvature at the 
crest is 1, the slopes of solutions A-C have similar shapes. 
In the traveling case (Fig. [5|, these rescaled slopes would ap- 
proach a limiting curve. But, for standing waves, oscillations 
develop, and a limiting curve does not emerge. 
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FIG. 5. Self-similar asymptotics of traveling water waves near the wave crest. (Left) Using a variant of our standing wave code, 
we computed 5 periodic traveling solutions with wavelength 2ir and particle speed q at the wave crest, chosen so that / = q 2 /2g 
has the values shown. We then plotted n(x)/l versus x/l, as well as the inner solution of [6|. The distance between successive 
crossings of the inner solution with its asymptote grows exponentially; thus, I must be extremely small to observe oscillatory 
behavior near the crest. (Right) As I — > 0, the slopes of the rescaled periodic waves approach the slope of the inner solution. 



It is instructive to compare our results to the traveling 
wave case. Longuet-Higgins and Fox [6J [7] showed that 
periodic traveling waves are asymptotically self-similar in 
two scaling regimes. If the wavelength, L, is held fixed 
as the crest tip sharpens, the limiting wave profile has a 
120° corner. This is the outer solution of [7], predicted by 
Stokes and proved to exist in [5] . If, instead, the fluid ve- 
locity at the crest remains fixed as the wavelength goes to 
infinity, the limiting wave profile is shown in Fig. [5j This 
inner solution crosses the asymptotes y = ±x/y/3 in- 
finitely often [B] , implying that traveling waves approach 
Stokes's limiting wave in an oscillatory manner, rather 
than monotonically, with L fixed. 

The oscillations in the standing wave case are of a com- 
pletely different nature. No choice of scaling will cause 
the curves in Fig. [4] to approach a limiting inner solu- 
tion. We believe these oscillations are caused by resonant 
modes in the two-point boundary value problem ([I]) with 
boundary conditions Q(x, ±T/4) = 0, treating T as a bi- 
furcation parameter. A resonant mode is a perturbation 
that nearly satisfies the linearized boundary value prob- 
lem. Such modes can be strongly excited in the process 
of computing standing waves, especially in finite depth 
[TBI |2"21 |2"3"] . Disconnections in the bifurcation diagram 
seem to occur when a resonant mode can be excited with 
more than one amplitude. For example, solutions I and 
J in Fig. [3] both contain a secondary, higher-frequency 
standing wave (the resonant mode) superimposed on a 
low-frequency carrier wave. The secondary wave sharp- 
ens the crest at J and flattens it at I, being 180 degrees 
out of phase from one branch to the other. 

We conclude that resonance is responsible for oscilla- 
tions and trumps self-similarity in determining the dy- 
namics of standing waves at small scales. This shows 
that, although under-resolved numerical simulations may 
exhibit self-similar dynamics, as happened in [15] . the 
true dynamics may be more complex. Recent work on 
singularity formation in free surface flow problems, such 
as droplet and bubble pinch-off [JJ |2] and wave breaking 
[3J, may also benefit from higher-resolution simulations, 
which could reveal new aspects of their dynamics. 
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